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Given a Taylor series with a finite radius of convergence, its Borel transform defines an entire 
function. A theorem of Polya relates the large distance behavior of the Borel transform in difi'erent 
directions to singularities of the original function. With the help of the new asymptotic interpolation 
method of van der Hoeven, we show that from the knowledge of a large number of Taylor coefficients 
we can identify precisely the location of such singularities, as well as their type when they are 
isolated. There is no risk of getting artefacts with this method, which also gives us access to some 
of the singularities beyond the convergence disk. The method can also be applied to Fourier series 
of analytic periodic functions and is here tested on various instances constructed from solutions to 
the Burgers equation. Large precision on scaling exponents (up to twenty accurate digits) can be 
achieved. 



U 



(N 
> 

in 
o 

O 
O 



> 

X 
J3 



INTRODUCTION 



In the late nineteenth Century, Pincherle [l[ and then 
Borel 0, [3| introduced what is now known as the Borel 
transformation: given a formal series in powers of the 
complex variable Z 



/(Z) = ^a„Z", 
one introduces the Borel transformed series 

F(c) ^ £ • 

n=0 

Since, for Re Z > 0, 

C"e-^« dC 



7.1 



it is useful to introduce the function 



^ ' ri=0 



(1) 



(2) 



(3) 



(4) 



which is formally the Laplace transform of F[C^ and 
which in this context is sometimes called the Borel- 
Laplace transform of F . 

Borel's motivation was predominantly to give a mean- 
ing to divergent series such as ^n!.Z" and the Borel 
transformation has been extensively used to resum diver- 
gent series appearing in physics (see e.g. Refs. 0, HQ)- 

In 1929 Polya observed that the Borel transforma- 
tion can also be used to obtain information about singu- 
larities of a Taylor series (in powers of 1 /Z) with a finite 
radius of convergence, in which case the function F{C) is 
entire. He proved a theorem relating the convex hull of 



singularities of f^^{Z) (the smallest convex set outside 
of which the function is analytic) to a function called 
the indicatrix of F{(), roughly the rate of exponential 
growth at infinity of F{(^) as a function of the direction 
(for precise definitions sec Section HV]) . 



Here we show that this theorem can be used in con- 
junction with high-accuracy numerical methods to ob- 
tain very precise information on singularities of Tay- 
lor and Fourier series. Singularities play an important 
role in fluid dynamics and condensed matter physics (see 
[liBESl and references therein). Using Polya's theorem 
to devise a practical numerical method would not have 
been possible without recent progress in high-precision 
numerical algorithms and, foremost, the new technique 
for asymptotic interpolation of van der Hoeven [ll] which 
can sometimes give remarkable precision (close to twenty 
digits) on scaling exponents. 



The paper is organized as follows. In Section [H] we re- 
call some known facts about Taylor and Fourier series and 
their singularities. In Section IlIII we give a presentation 
of van der Hoeven's method from an applied mathemat- 
ics point of view and show how it works in practice, using 
known results for the Burgers equation. In Section lTVl we 
give an elementary introduction to Polya's theorem. In 
Section|V]we present our new method, which we propose 
to call BPH (Borel-Polya-Hoeven), for determining the 
convex hull of singularities and, for the case of isolated 
singularities on this hull, their positions and type. In 
Section IVII we test BPH using again the Burgers equa- 
tion. In Section [VIII we discuss open problems and make 
concluding remarks. 
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II. FROM TAYLOR AND FOURIER 
COEFFICIENTS TO SINGULARITIES 

We first recall the close relation between Fourier and 
Taylor series for analytic functions. Let u{x) be a 2tt- 
periodic function which is analytic in some neighborhood 
of the real axis in which it can be extended to a function 
where z — x + iy. After subtraction of a suit- 
able constant we can assume that J^^ u{x)dx ~ 0. The 
Fourier-series representation of u reads 



u{x) = 



E 

k=±l, ±2, 



e ^''''u{x)dx 



(5) 



(6) 



We denote by u^{x) (resp. u~(x)) the partial sum of 
the Fourier series ([5]) with fc > (resp. fc < 0), which is 
analytic in the upper (resp. lower) half plane y > (resp. 
y < 0). Each of these two functions can be written as a 
Taylor series by an exponential change of variable: 



u+iz) 
u~ (z) 



k>0 

E 

A;>0 



UkZ\ 



u-kZ\ 



Z 



z 



(7) 
(8) 



This can be shown directly by applying standard steepest 
descent asymptotics to the integral © [11] • 

When the radius of convergence of a Taylor series is 
determined by a single singularity of the type assumed 
by Darboux, the knowledge of a sufhciently large num- 
ber of Taylor coefficients with enough accuracy permits 
an accurate determination of the position and type of the 
singularity. This can be done by an iterative algorithm 
developed by Hunter and Guerrieri or by the asymp- 
totic interpolation method discussed in Section [IIII 

Sometimes there are two Darboux-type singularities on 
the convergence circle or, equivalently, the periodic func- 
tion (z) has two singularities with the same imaginary 
part. The interference of the two singularities produces 
then a sinusoidal modulation of the Taylor coefficients. 
This can still be handled by an iterative algorithm 
but not directly by the asymptotic interpolation method, 
for reasons explained in Section [III Bl The BPH method 
of Section |V] can handle not only the case of two or 
more isolated singularities on the convergence circle but 
also "hidden" singularities located beyond this circle (or 
within this circle if the series is in inverse powers of Z) , 
whose contributions to the Taylor coefhcients are expo- 
nentially smaller than any term in (jlip . From an asymp- 
totic point of view these contribution are "beyond all 
orders" . 



Obtaining the singularities of an analytic periodic func- 
tion from its Fourier coefficients is just basically the same 
problem as obtaining the singularities of an analytic func- 
tion f{Z) = 'Yl^^o'^nZ^ from its Taylor coefficients a„. 
Hadamard's formula gives us the radius of convergence 
of the Taylor series, namely the distance to the origin of 
the nearest singularity (ies). If we happen to know that 
this is an isolated singularity at Z^,, we can relate the 
singular behavior near Zi, to the asymptotic behavior of 
the an by the Darboux theorem [13, [13, [i3|- For this 
one assumes that, in a neighborhood of Z*, the function 
f{Z) is given by 

.f(Z) = il-Z/Z,)-'r{Z)+a{Z), ^ 0, -1, -2 . .(9) 



riZ) = J2b,il-Z/Z.f 



(10) 



k=0 



where the functions r(Z) and a{Z) are analytic in some 
disk centered at the origin with a radius exceeding \Zi^\. 
It then follows that, for large n, 



fln-E 



k=0 



{-l)%kZ':~''T{n + iy^k) 
n\T{h' — k) 



(11) 



The leading term is simply a„ ~ hQ-n'^^^ / {Zl^)T{v). Ap- 
plied to the Fourier series ([5]), the leading-order Darboux 
formula can be recast as follows: a branch-point singular- 
ity with exponent —v of u{z) at a location in the lower 
complex plane implies that for k — > -|-oo the Fourier co- 
efficient Uk is asymptotically proportional to k'^~^e~^^^* . 



III. THE ASYMPTOTIC INTERPOLATION 
METHOD 

Suppose that we have a function G{r) of a scalar posi- 
tive variable r for which we suspect that it has, for large r, 
an asymptotic expansion with a leading term Cr~"e~^^ , 
as in the Darboux theorem (|lip . but that we only know 
its values numerically with high accuracy (tens to hun- 
dreds of known digits) on a regular grid tq, 2ro, . . . , NrQ 
with a large number N of points (from fifty to thousands, 
depending on the problem). We set 



Gn = G{nrQ), n=l,2,. 



.N 



(12) 



Can we determine parameters such as C, a and S with 
high accuracy? One way is of course just to ignore the 
subleading corrections and to try a least square fit of the 
data to the functional form Cr~"e~^'", after taking a log- 
arithm. One then has the awkward problem of having to 
pick a fitting interval of values of n; the procedure usually 
gives poor accuracy and the determination of subleading 
corrections is almost impossible. 

A better way, used for example in Refs. [H, [H, [i3| j is 
to notice that, if we take the second ratio, defined as 



p _ GnGn-2 _ 1 

''"=^rr"l (^^ 



(13) 



then both the constant C and the exponential drop out. 
Assuming then n to be sufficiently large that we can ig- 



3 



nore subleading corrections, we obtain 



In R„ 



a 



ln(l - 1)2) 



(14) 



The other two parameters C and S appearing in ([iD]) 
are then easily determined. If the remainder, that is the 
discrepancy between the value of a predicted by pi)) . 
which we denote a„, and its Umit aoo for n —f oo, tends 
to zero in a known functional way, e.g., exponentially or 
algebraically, then we can extrapolate the q;„'s to infinite 
values of n using, e.g., one of Wynn's algorithms [H, US] 
(see, Ref. 2l| for a review of extrapolation methods). 
We shall come back briefly to such issues in Scction llllBI 
Without knowing something about the functional form of 
the subleading corrections which control the remainder, 
extrapolation may not work very well because the choice 
of the appropriate algorithm depends on the functional 
form of the remainder. 

Recently, van der Hoeven introduced the asymptotic 
interpolation method ll| which allows in principle the 
determination of the asymptotic expansion of G„ beyond 
leading-order terms. When the function Gn is known 
with very high precision and up to sufficiently large val- 
ues of n, parameters such as the scaling exponent a can 
sometimes be determined with extreme accuracy, as we 
shall see in Section IlII Al An important feature of the 
asymptotic interpolation method is that it uses the de- 
termination of subleading terms to improve the accuracy 
on leading-order terms. 

Here we shall just give a short elementary introduc- 
tion to the asymptotic interpolation method for the case 
when the data G„ are real numbers. There are several 
variants of the asymptotic interpolation method; ours dif- 
fers occasionally from that of Ref. [11]. The basic idea of 
the asymptotic interpolation method is to perform sim- 
ple "down" transformations on the data G„ which suc- 
cessively strip off leading and subleading terms. After a 
number of such down steps which depends on the qual- 
ity of the data, the transformed data become sufflciently 
simple to allow a straightforward interpolation step. The 
list of down transformations which are needed is given 
hereafter. 



I Inverse: G„ — > 
R Ratio: G„ 



G„ 
G„ 



SR Second ratio: G,, 
D Difference: G„ — 



G„G„_2 
Gn ~ Gn-1 



At each stage, tests are applied to decide which of the 
four transformations should be applied in order to favor 
the stripping process as much as possible. If |G„| < 1 
for large n, apply I; otherwise proceed. If |G„| grows 
"slowly" at large n (we found that a useful operational 
definition is to see if the growth can be identified as al- 
gebraic with a rather well defined exponent), apply D; 



otherwise ("fast" growth), apply R. In addition, if |G„| 
grows or decreases exponentially at large n, we found 
that it saves time to apply SR; also, if |G„| is a slowly 
decreasing function, it is more convenient to apply — D. 
Note that the differences or ratios involved in R, SR and 
D are backward; this conveniently keeps the maximum 
index N fixed. 

When the procedure is iterated, after a while, an "in- 
terpolation stage" is reached where the data can be 
asymptotically interpolated in a simple fashion, typically 
by a constant plus a small remainder tending to zero at 
large n. Basically this means that we have successfully 
stripped off a certain number of terms in the asymptotic 
expansion. For the kind of data which we are considering 
here, the most useful interpolation stages usually arise 
at the sixth and thirteenth stages (counting the original 
data as stage zero). 

There are two effects which limit the number of stages 
which can be applied to a given set of data. First, when- 
ever a ratio or a difference are taken, the precision of 
the data (i.e. the relative rounding error) deteriorates; as 
the number of transformations applied increases, round- 
ing errors make the data increasingly noisy, beginning 
usually with the highest values of n. Second, the in- 
terpolation stages require sufficiently large values of N, 
since the constant asymptotic behavior at large n may be 
preceded by non-trivial transients. For a given resolution 
N and a given precision, the procedure must be stopped 
at the latest interpolation stage not significantly affected 
by the two effects just mentioned. In practice we should 
have a significant range of values of n over which the data 
are almost constant and not affected by rounding noise 
(if the rounding noise is very low this range may extend 
all the way to N). When the down process is stopped the 
data are interpolated and the process is reversed, by ap- 
plying "up" transformations which are the inverses of the 
down transformation in the reverse order. The inverses 
of the D, R and SR transformations involve one or two 
unknown additive or multiplicative constants which are 
determined using the highest known values of the G„ and 
of their down transforms. 

When the process is completed, the data are asymp- 
totically expressed as a truncated transseries. Roughly, 
a transseries is a formal asymptotic series involving in- 
teger or fractional powers, logarithms, exponentials and 
combinations thereof [l^, [l^, 111] . 

A worked example will now give the reader a more 
concrete feeling. 



A. Testing the asymptotic interpolation method on 
the Burgers equation with a single-mode initial 
condition 



Here and in Section IVII we shall perform tests using 
the one-dimensional inviscid Burgers equation 



dtu{t, x) + u{t, x)dxu{t, x) = , 



(15) 



4 



with a 27r-periodic real initial condition uq(x) having a fi- 
nite number of Fourier harmonics. We begin by recalling 
some well-known facts about the solution and the singu- 
larities of the inviscid Burgers equation. Eq. ([TS]) has an 
implicit solution in Lagrangian coordinates 



u(t, x) — UQ{a); X ~ a + tuQ^a) 



(16) 



Up to the time of the appearance of the first shock, 
the Lagrangian map a t—^ x has a Jacobian 



(17) 



J{t, a) = 1 + tdaUo{a) 



which does not vanish in the real space domain and ([T| 
defines a unique real solution. This solution has singular- 
ities in the complex domain (with the real part defined 
modulo 2tt) at locations which are the images by the 
Lagrangian map of the zeros of the Jacobian J. Generi- 
cally these are simple zeros. The singularities in Eulerian 
coordinates are then square-root branch points. The so- 
lution can also be written explicitly using the Fourier- 
Lagrangian representation [29l . Isoj , which in a special 
case was actually discovered earlier by Platzman [3l| . In 
the periodic case, the simplest representation, called the 
third Fourier-Lagrangian representation, valid for k ^ 0, 
is 



u{t,x) = J2 ^'""^kit) 



(18) 



k— — c 



Uk{t) = - — 



1 



2i7rfci 



2tv 



e-iHa+tuo(a))^^ ^ (19) 



In this section we take the "single-mode" initial condi- 
tion 



(20) 



Mo(a) = - -sina , 



for which the first real singularity is at t^, = 2. Using Jit 
and the integral representation of the Bessel function J„ 
of integer order n (see, e.g., Ref. [s^l p. 360), one finds 



(21) 



For convenience, we shall consider the solution att = 1. 
This single-mode solution has only one pair of complex 
conjugate singularities on the imaginary axis at 



±i(5, 



6 = In 2 



(22) 



Bessel functions of large order and arguments have 
an asymptotic expansion (in the sense of Poincare) , ob- 
tained through the method of steepest descent by De- 
bye I33l. (The matter is also discussed in Chap. VIII of 
Ref. [34|). Debye identified various asymptotic regimes 
which, in our notation, depend on whether t is less or 
larger than t^, and is or is not very close to i*; his classi- 
fication is in one-to-one correspondence with that of the 



various regimes relating to preshocks, as discussed, for 
example, in Refs. [1^, [s^]. When t — 1, well before t^,, 
the relevant Debye expansion for k — > +oo is: 



-.k- 



\W5 

where S is given by 
3$ - 



E 



Jn{2/V3) 



A:" 



(23) 



71 (0 

72(0 
73(0 



24 ' 
81^2 - 462^4 385^6 

1152 ' 
30375^3 _ 369603^5 ^ 765755^^ - 425425^^ 



414720 



(24) 



and the higher-order polynomials 7n(C) satisfy recurrence 
relations given, e.g., in Ref. [l^l- The leading term of this 
expansion follows also from the Darboux theorem (fTTjl . 

Let us now show that the asymptotic interpolation 
method, as outlined in Section IIIIl when applied to the 
Fourier coefhcients of the single-mode solution (PT|) can 
recover a suitably truncated version of the Debye ex- 
pansion We use all the Fourier coefficients with 
k — 1,...,N, where N ~ 1000 and define our initial 
data set as G„ = -u„(l)/i. 

Each coefficient is calculated with an 80-digit preci- 
sion (using Mathematica® and 120-digit working pre- 
cision). The basic transformations and their in- 
verses are implemented numerically in 80-digit pre- 
cision, using the high-precision packages GMP and 
MPFR available from http : / / www . swox . com/ gmp/ and 
http : //www . mpf r . org/ . 

With these data we are able to reach stage 13. The list 
of successively applied transformations, resulting from 
the tests given in Section Hill is 



SR, -D, I, D, D, D, D, I, D, D, D, D, D (25) 

Fig. [T] shows the first six stages. It is mostly intended to 
bring out overall features and to make clear which of the 
four transformations is to be selected at the next stage 

It is very easy to understand why the first six stages are 
as listed above. Indeed, let us suppose that, to leading 
order, G„ = Cn-^e""^". We can work out analytically 
the various transforms and we list hereafter the result up 
to the sixth stage, displaying only the leading and when 
needed the first subleading term in the large-n expansion: 



Cr 





^ u 


a 


-D 










D 3n 


D 


3 


'2^ 


a 




a 



2a 

71 



2^ 



(26) 



It is seen that stages 1 and 6 are interpolation stages at 
which the data are asymptotically fiat. Stage 6 is par- 
ticularly important since the asymptotic value 3/a gives 
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FIG. 1: Numerical output from asymptotic interpolation at stages 1-6. Stages 1-5 are represented in linear coordinates. For 
stage 6 we represent the difference between Ge{n) and its asymptotic value 2 in lin-log coordinates. 



the scaling exponent a. According to ([23|) . for the Burg- 
ers single-mode solution, the asymptotic value should be 
2. 

Let us now show in some detail how the asymptotic 
interpolation technique works to give us the asymptotic 
expansion of G„. We begin by limiting ourselves to a 
six-stage procedure. The successively transformed data 
wiU be denoted G^^^ . . . , G'-^\ Following Ref. [ll|, we 
interpolate G*^^-* by 3 /a. How cleanly this can be done 
is visible in the last of the graphs in Fig. [1] where we 
show the discrepancy between gI^'' and its asymptotic 
value 2. This discrepancy falls to about 10^^" at the 
upper end of the range. Then we determine G^^-* by 
inverting the relation G'^' = DG'^^-'. This involves an 
unknown additive constant which is determined from the 

(5) 

last data point G)^ . Then we continue inverting the 
D operators appearing at stages 4 and 5, each time us- 
ing the last point to obtain the additive constant. In 
this way we obtain a cubic polynomial for G*^'^^ We 
then invert the operator I and obtain the inverse of the 
aforementioned cubic polynomial, which can be written 
—2an~^{l + di/n + d2/n'^ + d3/n^ + . . .) with in principle 
well defined constants di, ^2, etc. Then we invert the op- 



erator SR; this can be done by taking a logarithm which 
will transform second ratios into second increments. At 
the end of the process we obtain the asymptotic expan- 
sion 



G„ 



Cn 

71 



72 , 73 
^ ^ + ~3 



O 



(27) 



It is actually simpler to start with (|27p and to apply 
successively the first six transformations listed in ([^5]) to 
identify the parameters. With the six-stage procedure 
we obtain G, a, 5, 71, 72, 73. Their values are given in 
Table HI 

It is seen that the coefficients G, a and 5 appearing in 
the leading term have a precision of at least 10~^". The 
precision of the coefficients 7^ for the subleading terms 
deteriorates with the order. 

We now turn to the analysis using a 13-stage proce- 
dure. This allows a much more precise determination 
of the aforementionned coefficients and, in principle, the 
determination of six additional terms in the expansion 
(P7)l . After stage 6, the next interpolation stage is stage 
13. This is easily shown by observing that the discrep- 



6 





a 


(5 


C 


6 stages 


1.49999999993 


0.4509324931404 


0.4286913791 


13 stages 


1.49999999999999995 


0.450932493140378061868 


0.4286913790524959 


Theor. value 


3/2 


0.450932493140378061861 


0.42869137905249585643 




71 


72 


73 


6 stages 


—0.17641252 


0.17295 


—0.401 


13 stages 


—0.17641258225238 


0.172968106990 


—0.406446182 


Theor. value 


—0.176412582252385 


0.1729681069958 


— 0.40d44d1o02 




74 


75 


76 


13 stages 


1.384160933 


-6.192505762 


34.5269751 


Theor. value 


1.3841609326 


-6.1925057618568063655 


34.526975286449930956 



TABLE I: Solution to the Burgers equation with single-mode initial condition: comparison of theoretical values with 6-stage 
and 13-stage asymptotic interpolation values for the first six coefficients in Debye's solution (|23|l . For 13-stage asymptotic 
interpolation we also give some of the higher order coefficients. 



ancy between Gn and its asymptotic value 3/a — 2 
is 0{l/n*), because all the lower-order terms have been 
stripped off by the first six transformations. More specif- 
ically, we have 

nP n" n' nP n-' 
C = O (^^) . (29) 

Stages 7-13 gives us the coefficients ci,..., cg and al- 
low us to find the remainder r„, as defined in We 
found that r„, determined by this 13-stage procedure, 
falls to about 5 x 10~^^ at the end of the range. As 
shown in Table HI the precision on the first six coefficients 
in the asymptotic expansion has improved very much and 
is now of a few lO^^'' for the exponent a. 

B. Further remarks on asymptotic interpolation 

The method of asymptotic interpolation is still in the 
development stage; improvements and new features are 
thus to be expected. Some are already suggested in the 
initial publication [llj. One rather straightforward ex- 
tension is from real to complex data. For rapidly grow- 
ing data, one can use logarithms instead of ratios. In 
Ref. [m it is recommended to take ratios or logarithms 
as often as possible and to define "slow growth" as slower 
than, say, n^/"^. This helps in identifying the functional 
form of the transseries expansion. Once this is known, we 
found that the values of the coefficients can be generated 
more efficiently by using a rather broad definition of slow 
growth, namely well-identifiable polynomial behavior. 

A very important issue is to determine how many 
stages are feasible with a given resolution TV and a given 
precision. We have found that the successive interpola- 
tion stages, at which the data are asymptotically fiat, 
have this flat regime preceded by longer and longer tran- 
sients. To make this more concrete we have investigated 



how far it is necessary to go to be within five per cent of 
the asymptotic value for the Burgers single-mode prob- 
lem. For stage 6 the asymptotic value is 2 and the data 
are within five per cent everywhere. For stage 13, the 
asymptotic value is about 0.33836513 and less than five 
per cent discrepancy holds for n > 25. The next inter- 
polation stage has number 20. The asymptotic value is 
2/7 but less than five per discrepancy holds only beyond 
n = 1620. In Fig.[2]we have represented the data at stage 
20. It is seen that the discrepancy is enormous until we 
reach well beyond n = 1000. Obviously, if N is not large 
enough the stripping of subleading terms performed by 
the successive stages must be stopped however high the 
precision of the data may be. A related issue is discussed 
at the end of Section 7 of Ref. 11|. 

Since rounding errors increase with the stage number, a 
certain balance must be kept between resolution and pre- 
cision. To investigate this quantitatively on the Burgers 
single-mode problem, we have artificially degraded our 
precision by adding random noise of various strengths. 
It appears that we need at least 16 significant digits at 
stage 6 and 27-35 significant digits at stage 13. 

If we are only interested in obtaining an accurate de- 
termination of a few terms in the expansion (|27[) , we may 
be able to retrieve them using asymptotic interpolation 
stopped at the sixth stage and continuing with a different 
strategy. Indeed we observe that, at the sixth stage, the 
data given by (p9| have the form 

g1^) = s + r„ , (30) 

where the remainder r„ decays to zero as n — > oo. This 
is a well studied situation in the theory of convergence 
acceleration by sequence transformations, whose goal is 
to replace the sequence ([5(11) by a transformed sequence 
having the same limit s but a much faster decaying re- 
mainder (see, e.g., Refs. [2l|,[2^). A simple and very pop- 
ular acceleration method appropriate for (P^)) is Wynn's 
rho-algorithm [l^ , although more sophisticated methods 
are known [2ll|. In our case it gives the correct value 2 
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FIG. 2: Interpolation stage 20 has the asymptotic value 2/7 as shown in the inset. The main figure shows the discrepancy 
G2o(n) - 2/7. 



with a 20-digits precision. This is even better than the 
13-stage asymptotic interpolation. Note that the choice 
of a particular convergence acceleration method depends 
crucially on the functional form of the remainder. With 
asymptotic interpolation this form can be determined 
rather than having to be assumed. Here a caveat is in or- 
der: if the data are not sufficiently asymptotic the mixed 
procedure just described will not work, for example be- 
cause the remainder has not yet settled down to algebraic 
decrease. A situation of this type seems to be present in 
the work on short-time asymptotics discussed in Ref. : 
the rho-algorithm does not improve the quality of the 
scaling exponent controlling the divergence of the vortic- 
ity at the singular manifold and much higher resolutions 
are probably needed for that problem. 

E.J. Weniger (private communication) has pointed out 
that asymptotic interpolation and sequence transforma- 
tions have technical features in common. In asymptotic 
interpolation one tries to annihilate leading terms in the 
asymptotic expansion, whereas in sequence transforma- 
tion one tries to shrink the remainder by annihilating 
its largest contributions, but the transformations used 
in both instance are often the same, for example, finite 
difference operators. 

A powerful method of asymptotic series analysis, 
widely applied in statistical physics, is the method of 
differential approximants, which can be viewed as a 
generalization of the Dlog Pade method 0]. In par- 
ticular it has been used to analyze self-avoiding walks 
(SAWs) and polygons (SAPs). We have applied the 
asymptotic interpolation method to data available at 



http : //www.ms . unimelb . edu. au/~iwan/polygons/Pol 
ygons_ser.html. The goal was to see how well we can 
reproduce the asymptotics of the number of self-avoiding 
polygons with 2n steps on square and honeycomb lattices 
psL . When analyzing the square lattice data for 

the largest available range, that is n up to 55, we found 
that asymptotic interpolation gives the value of the criti- 
cal point correct to 9 decimal places, whereas differential 
approximants give about 3 additional digits. We observe 
that (i) the actual implementation with asymptotic in- 
terpolation is somewhat simpler and (ii) asymptotic in- 
terpolation is not limited to problems which can be well 
approximated by solutions of low-order linear differential 
equations. 



The method of asymptotic interpolation is, in our opin- 
ion, very useful but is of course not the panacea. One dis- 
ease it cannot directly cure is the presence of sinusoidal 
oscillations. For example if the analytic function f{Z) 
has two complex conjugate singularities at i?e*"^* on its 
circle of convergence, large-order Taylor coefficients will 
present a sinusoidal oscillation with a wavelength propor- 
tional to (t>i,. After any number of stages, this oscillation 
is still present and the data cannot be interpolated by 
a constant. As we shall see now, a Borel transforma- 
tion takes care of this problem and can also bring hidden 
singularities to the foreground. 
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IV. POLYA'S THEOREM 

Here we just want to give the reader a good feeling of 
what the theorem states and a heuristic derivation. We 
begin with examples discussed in Section 32 of Ref. ■ 

Let c = |c|e~''*' be a complex number and consider the 
function 



F(C)=e^^ = l + ^ + ^ + ... , 

which corresponds to the choice a„ = c" in ([2|). 
Borel-Laplace transform, given by (0]), is 



1 



1 



Z-c 



(31) 
The 

(32) 



It has a pole at Z — c, whereas -F(C) is an entire function 
(analytic in the whole complex domain). We set ( = 
re"^ and let r oo, holding the direction fixed; the 
modulus of F(C) — el'^l'"'''* ^, in the direction 0, varies 
exponentially at the rate h{(j)) = |c|cos(0 — 7), called 
the indicatrix of F. We define k{(j)) = Re (ce~"^) = 
|c| cos((/) + 7) . This is the (signed) distance of the origin 
to the line normal to the direction passing through the 
pole c and is called the supporting function of the (single) 
singularity. We observe that 



H^) = ki- 



rn 



This relation is the simplest instance of Polya's theorem. 

Next, following again Polya's Section 32, we want 
to have n distinct poles at the complex locations 




Cl, C2, 



For this we take complex linear combi- 3' Construction of the supporting function fc(0) of the 



nations with non- vanishing coefficients Ci, C2 



Cri 



F(C) - Cie^i^ + C2C^'^ + ... + Cpe^^^ . (34) 
The Borcl-Laplace transform is 
Cl C2 



f^{Z) = 



Z — Cl Z — Cl 



Cp 

Z — Co 



(35) 



For any G [0, 27r], we define now the indicatrix and the 
support function, a little more formally, as 



h{^) EE limsupr-Mn|F(re''^) | , 



(36) 



(in the present example, the limsup is just an ordinary 
limit) and 



A:(0) = sup (X cos ( 



ysin0) = sup {Re (Ze"''^)} , 

(37) 

where Z ~ X and K = {ci, C2, . . . , Cp} is the sin- 

gular set. Since there is a finite number of singularities, 
the sup operation is just the same as the maximum. We 
define a supporting line of as a line which has at least 
one point in common with K and such that all the points 
of K are in the same half space with respect to the line. 



set K of singularities of / {Z). The singularity c\ is on the 
convergence circle; the convex hull of K is defined by Ci, C2 
and C3. The singularity C4 is inside the convex hull. 



The intersection of all these half spaces is the convex hull 
of K . In the present case this is just the smallest convex 
polygon containing all the poles. It is readily seen that 
A:(0) is the (signed) distance of the origin to the support- 
ing line normal to the direction (see Fig. |3|). The rate 
of growth of F(C) in the direction is obviously that of 
the fastest growing of the p exponentials in p4|) , which 
is precisely k{—(l)), so that Polya's relation ([55)1 holds. 

He proved a much more general theorem: Let f^^{Z) 
be an analytic function defined by the Taylor series ([4]) 
in powers of l/Z which has a finite non-vanishing ra- 
dius of convergence H , and let K be the smallest convex 
compact set containing its singularities, then (i) the Borel 
transformed series Q defines an entire function of expo- 
nential type, and (ii) the indicatrix h{(f>) of F{C,), defined 
by p6p . and the supporting function k((j)) of K , defined 
by (|37p . are related by h((j)) = k{—<j)) and H — sup^ h{(f)). 

(An entire function F{Q is said to be of exponential 
type if its modulus is bounded by Ae^'^l, where A and a 
are suitable positive constants.) 
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Polya's proof (not given here) makes use of the fact 
that f^^ and F are Laplace transformed of each other, 
specifically, 

poo 

/^i^(Z)= / F{Oe-^<dC, (38) 

-| /•a+i oo 

F(C) = — / /Bi^(Z)e^«dZ, (39) 

where a is any real number such that the singular set K 
is entirely contained in Re Z < a. 

Observe that no particular assumption is made regard- 
ing the type of the singularities which can be isolated (e.g. 
poles or branch points) or continuously distributed (nat- 
ural boundary). Inside the circle of convergence \Z\ = H, 
the series (|3]) is divergent. However if the whole circle is 
not a natural boundary, the function f^^{Z) can be an- 
alytically continued to suitable Z's inside this circle and 
the pair of integrals can be viewed as a way of 

resumming the divergent series (U). 

In applications it frequently happens that all the "edge 
singularities" , that is, those determining the border of the 
convex set K are isolated. This border is then piecewise 
linear, as in the case of n poles discussed above. The 
angular dependence of the supporting function is then 
given by k{(j>) = \cj\ cos((/) + 7^) in the angular interval 
< (f) < 4>j for which the supporting line normal to 
(j) touches K at cj — \cj\e~^^^ (see Fig. [3]). If k{(j)) is 
known with high accuracy, then the positions of the edge 
singularities CjS can also be determined accurately. 

Moreover we can then determine the type of an isolated 
singularity at Cj by studying the asymptotic behavior 
of F(C) along rays ( = re^'^ with large r, in a suitable 
angular interval. For example, let us assume that, near 
Cj the function f{z) has an algebraic singularity and is 
to leading order proportional to {Z — Cj)"^^, where the 
exponent a is real and not a positive integer. (If a — 
1 > 0, this behavior is assumed for a suitable first- or 
higher-order increment of /.) After shifting the contour 
of integration to follow the boundary of K near cj (cf. 
Fig. S]), apphcation of steepest descent [131 to ([M)) with 
C = re"^ taken in the angular sector (jij-i < ~4' < 0i 
and r — > 00 yields 

G(r) = |i^(re''^)| = Cr^°e''('^)'-[1 + £(r)], (40) 
M0) = |c,|cos(0-7,)- (41) 

Here C is a positive constant and e(r) tends to zero for 
r ^ cxD at a rate which depends on what is assumed for 
subleading corrections to the (z — Cj)"~^ singular behav- 
ior. If we are able to identify the algebraic prefactor to 
the exponential in (|40p. we can find the exponent a of 
the algebraic singularity. 

Non-algebraic singularities can be handled similarly. 
For example, if near Cj the function f{z) behaves as 
Q^/iz-cj) ^ application of steepest descent shows that in- 
stead of the algebraic prefactor proportional to r~" which 
appears in (j40p . we obtain an exponential prefactor pro- 



a + ioo 



n 




a - ioo 

FIG. 4: Contour of integration for computing the inverse 
Laplace transform of f^^{Z) (dashed line) and its deforma- 
tion to obtain the asymptotic contribution from the singular- 
ity at Cl (continuous line with arrows). 

portional to e^'^<^°'^i'P/'^)Vr _ Furthermore, if all the singu- 
larities on the convex hull of K are isolated, then ([77)) 
remains valid: the indicatrix is piecewise a cosine func- 
tion. 

How this is done in practice will be discussed in the 
next section. 

V. THE BOREL-POLYA-HOEVEN METHOD 

As we have seen in Section Hill the asymptotic inter- 
polation method, applied to the Taylor coefficient of an 
analytic function with a finite radius of convergence de- 
termined by a single isolated singularity allows one to 
identify its location and type. This is not the case if 
there is more than one singularity on the convergence cir- 
cle. Furthermore "hidden" singularities are not directly 
retrievable from the asymptotics of Taylor coefficients. 

We can however take advantage of Polya's theorem 
f Section livp to replace the analysis of large-order Taylor 
coefficients of an analytic function by the analysis of the 
behavior of its Borel transform at large distances in the 
complex C plane along various rays. This behavior can 
be found by asymptotic interpolation, from which we can 
then construct the convex hull K of the singularities and 



10 



obtain their type when they are isolated. This is the 
BPH strategy which we now describe in a more detailed 
way. 

We start from a truncated Taylor series in inverse pow- 
ers of Z 

N 

/tH^) = E^' (42) 

71=0 

with N terms, each of the coefficients being known with a 
precision e. We construct the associated truncated Borel 
series 

N 

^t(C)-E^C"- (43) 

71=0 

We choose a certain number of discrete angular di- 
rections characterized by their angle 0. Along each ray 
C = re"^, we evaluate the modulus of the Borel series at 
M points spaced by a constant distance (mesh) tq: 

G™ = |FT(mroe'^)|, m = 1, 2, . . . , Af . (44) 

Then we apply the asymptotic interpolation method of 
Section [nT] to identify a large-m leading-order behavior. 
For example, for algebraic singularities we have 

Gm ^ C(</.)(mro)^"('^'e''('^'""^° . (45) 

This gives us the constant C{4>), the prefactor exponent 
— a((/)) and the indicatrix h{(p) for the discrete set of di- 
rections. Polya's theorem then gives us the supporting 
function A:((/)) — h{—(t)) of the set K of singularities of 
the Taylor series. As we have seen in Section IIVI if the 
singularities on the convex hull of K are isolated and 
are located at |cj|e~''''^ , then the supporting function is 
piecewise a cosine function, given by \cj \ cos{(j) -\- ^ j) . The 
exponent a gives us the type of the singularity: a branch 
point (or a pole) of exponent a — 1. Other types of sin- 
gularities, for example of the exponential type discussed 
near the end of Section |Vl are handled similarly after 
identification of the appropriate asymptotic behavior. 

In practice, we have to choose the set of discrete direc- 
tions, the mesh tq and the maximum number of points 
M on each ray. If we happen to know the number p of 
isolated singularities and, at least approximately, their 
positions we can pinpoint the latter by taking 2p suit- 
able (j) directions. This is however rarely the case. We 
recommend taking a fairly large set of directions (for ex- 
ample 500 uniformly spaced directions) in order to reduce 
the risk of missing one or several of the cosine functions. 
The natural choice for the mesh rp is where H is the 
radius of convergence of the Taylor series. An approxi- 
mate value is i?approx = (I/"-) In |a„| for large n, which 
is roughly constant. For the determination of the largest 
distance rmax = Mtq we limit ourselves to the case where 
the function F(re^'^) grows at large distances, that is 
h{(j)) > (otherwise there are severe numerical prob- 
lems). r,„ax is then determined by the condition that the 



last term unC^/NI of the (truncated) Borel series (|43p 
should introduce a relative error in the determination of 
F(re' '^) which does not exceed the precision e with which 
the Taylor coefficients are known. A rough estimate for 
IflAfl is and for |F(rniaxe''^) is e'''"--" . Using the Stir- 
ling formula, we find that, to leading order M ~ N (the 
dependence on e appears only in subleading logarithmic 
corrections) . 

We mention that an improvement would be to replace 
a mere polynomial truncation of the Borel series by a 
suitable resummation/acceleration method for comput- 
ing entire functions [38] . This could be crucial for deter- 
mining negative indicatrix values, that is, when i^(C) is 
exponentially decreasing at large C- 

It is of interest to know how well we can separate two 
discrete singularities. By Polya's theory, each singularity 
contributes an exponential term to the large-r behavior 
of the modulus of the Borel transform. If r becomes 
sufficiently large compared to the difference in the two 
e- folding rates, only one of the two singularities will be 
seen. By suitably changing the direction of the ray in the 
C-plane we can then focus separately on each singularity. 
The worst case for discrimination is when we have two 
singularities which are at the same distance of the ori- 
gin. Assuming that this common distance is comparable 
to the radius of convergence H and denoting by A the 
distance of the two singularities, we find that the largest 
discrepancy in e-folding rate is roughly A^/if. Denot- 
ing, as above, by M the maximum number of point on 
a ray, we find that good separation requires the separa- 
tion parameter AIA^/H^ to be large. Since discrepancies 
are amplified exponentially, a separation parameter of 10 
may suffice. 

We shall not here discuss issues of algorithmic com- 
plexity, such as the reduction of the number of operations 
to evaluate the truncated Borel series. In applications 
the complexity of the numerical calculations needed to 
accurately determine the Taylor coefficients will usually 
exceed very much what is needed for the BPH analysis. 



VI. TESTING BPH ON THE BURGERS 
EQUATION WITH MULTIMODE INITIAL 
CONDITIONS 

To test the BPH method we need a Taylor series having 
either a pair of singularities on the convergence circle or 
"hidden singularities" . As in Section IIII Al this can be 
done using 27r-periodic solutions of the inviscid Burgers 
equation. The 2-mode initial condition 

uo{a) = Al sin a -I- A2 sin(2a) , (46) 
Ai = -l/2, A2 = (l/16)(4- v/(14)) + e, 
e = 1/150 
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produces at t = 1 a solution m(1, z) having, in Eulerian 
coordinates, singularities at 

z± = ± 0.1103542160016972443 

± i 0.737097018253664793 . (47) 

Henceforth we shall concentrate on the singularities of 
u{l, z) in the lower half plane, which are also the singular- 
ities of the function u~^{l, z), the sum of Fourier harmon- 
ics with fc > 0. Note that there are two singularities with 
the same imaginary part and opposite real parts (this is a 
consequence of the symmetry a t—^ —a, uq i— s- —uq of the 
initial condition and of the complex conjugate symme- 
try). When the Fourier series for u+(l, z) is transformed 
into a Taylor series in inverse powers of Z by setting 
Z = c^'^, the z singularities get mapped onto two com- 
plex conjugate Z singularities 

Z± = 0.4755903313336372343±i 0.0526974896343733942 . 

(48) 

The 3-mode initial condition 

uo{x) — Ai sin a; + A2 sin 2x + A3 sin 3a; (49) 

Ar = A A.^i^ + 1, 
^2' 16 50' 



A. = 



1 

Too 



(50) 



produces at t = 1 a solution u{l,z) having, in Eulerian 
coordinates, in the lower complex half plane singularities 
at 

zi* = -i 0.4608974136239120258 (51) 
zf^ = ±0.8575677577466957833 

-i 1.1175132271503113898. (52) 

The zii, singularity is on the imaginary axis and is the 
closest to the real domain. The other two are further 
away (hidden) and symmetrically located with respect 
to the imaginary axis. The corresponding Z singularities 
are 



Zi. 



^2* 



: 0.6307173770893952917, (53) 
: 0.2140094820693456182 
±i 0.2473645913888956747 . (54) 



are shown in Fig. O 

To apply the BPH method we generate the Fourier 
harmonics with k = 1, . . . , N = 1000 using the third 
Fourier-Lagrangian representation ([TO)) . The Lagrangian 
integrals are again calculated with 80-digit precision and 
120-digit working precision. 

The Borel transform is calculated for 20 values of 
between and tt/2 (a symmetry (j) — > —0 makes it un- 
necessary to take (p < and for > 7r/2 the indicatrix 
is negative). For the mesh we take tq = 1. The total 
number of points on each ray is M = 500. Along each 
selected ray, application of six-stage asymptotic interpo- 
lation gives us C, a and h. We can check that the indi- 
catrix is piecewise a cosine function as implied by ([37]). 
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FIG. 5: Positions in Z-p\ane of the singularities for the 3- 
mode initial condition. Continuous line: circle of convergence; 
dashed line: image of the real domain by the map 2 e~'^. 



Least square fits allow us to identify the parameters of 
these cosine functions and thereby to find the locations 
of the singularities. We recover the known values with an 
accuracy of about 10~^. We found that the accuracy on 
"hidden singularities" is comparable to that on directly 
visible ones. We also found that N — 1000 is not suf- 
ficiently asymptotic for a 13-stage analysis of the kind 
described in Section UlIAl 



VII. CONCLUDING REMARKS 

One central theme of this paper is the use of a Borel 
transform, in conjunction with Polya's theorem, to re- 
veal singularities not directly accessible from the asymp- 
totic behavior of the Taylor/Fourier coefficients. A very 
useful property of the Borel transform of a Taylor se- 
ries (in inverse powers of Z) is that its large-distance 
behavior encodes information not only about those sin- 
gularities of the Taylor series located on its convergence 
circle, but also about other singularities "hidden" inside 
this circle. Actually the Borel transform, followed by a 
Borel-Laplace transformation is a way of performing an- 
alytic continuation. Recovering hidden singularities from 
a Taylor series has important applications in a number of 
fields; many of the known techniques have been reviewed 
by Guttmann f^. 

To the best of our knowledge Polya's theorem has never 
been used as a numerical tool for identifying singulari- 
ties. The theorem is of a very general nature and as- 
sumes nothing about the nature of the singularities; this 
has the great advantage that we do not have to distin- 
guish between true and spurious singularities, as is the 
case, for example, when using Pade approximants and re- 
lated methods. The principal drawbacks are that (i) not 



12 



all hidden singularities are accessible, only those located 
on the convex hull of the singular set, (ii) pairs or clus- 
ters of singularities situated to close to each other may 
not be easily distinguishable, and (iii) enough terms in 
the series must be known to be able to actually obtain 
the asymptotic behavior of the Borel transform. When 
hundreds to thousands of Taylor coefficients are known, 
alternative mathematically well-founded techniques may 
become competitive, for example the old Weierstrass an- 
alytic continuation method; thanks to recent algorithmic 
discoveries it can be performed quite efficiently 39, 40]. 

The other theme of this paper is the asymptotic in- 
terpolation method of van der Hoeven which is here 
used both directly (when Darboux's theorem is applica- 
ble) and indirectly by means of Polya's theorem. When 
a large number of Taylor/Fourier coefficients are know 
with sufficient accuracy, asymptotic interpolation can 
give truly remarkable results, providing us not only with 
very accurate leading terms but also with several sub- 
leading corrections. As we have seen in Section Hill there 
is usually a well-defined relation between the number of 
subleading correction terms and the number of stages 
of the procedure which can be achieved. The latter de- 
pends crucially on the number of known coefficients and 
on their precision. For example if the data have only 
double precision, it is unlikely that more than six stages 
can be achieved. Asymptotic interpolation might than be 
viewed as an overkill compared to more standard tech- 
niques, but it is worth stressing that asymptotic interpo- 
lation is very easy to implement. 

Which kinds of problems are most likely to fall 
within the prongs of full-strength Borel-Polya-Hoeven- 
type analysis? This depends crucially on the computa- 
tional complexity of the problem, that is the dependence 
of CPU requirement and storage on the number of co- 
efficients N. As pointed out by Guttmann phase 
transition problems formulated on a lattice require usu- 
ally enumerating diagrams and the number of these tends 
to grow exponentially with order, while fluid dynamics 
problems generally have only polynomial complexity. In 
connection with phase transitions our BPH method is 
likely to be less precise than alternative methods such as 
differential approximants, but it can usefully supplement 
them to ascertain that the singularities identified are not 
artefacts. 

In fluid dynamics one outstanding problem is the issue 
of finite-time blow-up for the three-dimensional incom- 
pressible Euler flow with smooth initial data [l^, [ll| . For 
initial data having simple trigonometric polynomial form, 
one can determine numerically a number of coefficients of 
the Taylor expansion in time of the enstrophy (integral 
of one- half the squared vorticity). This was done for the 
Taylor-Green flow by Brachet et al. [l^l (yielding 40 non- 
vanishing coefficients calculated with quadruple working 
precision) and for the Kida-Pelz flow by Pelz and Gu- 
lak [4^ (yielding 16 non-vanishing coefficients having at 
least 40-digit precision). Because the number of coeffi- 
cients is rather small, there is no consensus on what the 



results imply for blow-up. Such calculations have a com- 
plexity 0{N^) which can however be reduced to 0{N'^) 
(up to logs) using the method of relaxed multiplication 
[39|. It is likely that a state-of-the-art calculation for 
flows with simple trigonometric polynomial initial condi- 
tions can give up to several hundred non- vanishing Tay- 
lor coefficients of the enstrophy, with a working preci- 
sion of several hundred digits. Another problem which 
can be tackled by series analysis is the analytic struc- 
ture of the two-dimensional incompressible vortex sheet 
(Kelvin-Helmholtz instability). It is known that an ini- 
tially analytic interface will develop a singularity in its 
shape after a finite time. Moore has made a prediction 
regarding this singularity which has been studied by 
various numerical techniques [45j . Again there is no 
consensus on the type of this singularity. 

It is of course of interest to extend to several di- 
mensions the BPH method, here presented only in the 
one-dimensional case. We observe that there exist 
multi-dimensional generalizations of the Borel transform 
[HiE^I and that the the asymptotic interpolation method 
can also in principle be extended to several dimensions 
[ll[. In several dimensions, singularities are not point- 
like; they reside on extended objects such as analytic 
manifolds and can have a much more involved structure 
than in one dimension. It is possible to partially recon- 
struct such objects using a variant of BPH. Furthermore, 
we note that Polya's theorem has been extended to sev- 
eral complex dimensions 48] (it is then referred to as the 
Ivanov-Stavskii theorem). Information on singularities 
can then in principle be obtained numerically in a way 
analogous to what has been done in Section |Vl This will 
be discussed elsewhere. 
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